Probing the nuclear EOS with fragment production 
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■ We discuss fragmentation mechanisms and isospin transport occur- 

, ring in central collisions between neutron rich systems at Fermi en- 

CSJ I ergies. In particular, isospin effects are analyzed looking at the cor- 
relations between fragment isotopic content and kinematical prop- 

4^ erties. Simulations are based on an approximate solution of the 

i | Boltzmann-Langevin (BL) equation. An attempt to solve the com- 

O ' plete BL equation, by introducing full fluctuations in phase space 

2 ■ is also discussed. 
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1 Introduction 

o ■ 

> ' 

In the last few years the increased accuracy of the experimental techniques 
rS [ has renewed interest in nuclear reactions at Fermi energies. Exclusive measure- 

ments, event-by-event analysis, and a 47r coverage allow a deeper investigation 
of the evolution of the reaction mechanisms with beam energy and centrality. 
New insights into the understanding of the nuclear matter equation-of-state 
(EOS) were gained [1] . In particular, recent experimental and theoretical anal- 
yses were devoted to the study of the properties and effects of the symmetry 
term of the EOS (iso-EOS) away from saturation conditions [2,3]. In partic- 
ular, heavy ion reactions with exotic nuclei at Fermi energies can be used to 
study the properties of the symmetry term at densities below and around the 
normal value. In central collisions at 30-50 MeV/A, where the full disassembly 
of the system into many fragments is observed, one can study specifically prop- 
erties of liquid-gas phase transitions occurring in asymmetric matter [3,4,5,6]. 
For instance, in neutron-rich matter, phase co-existence leads to a different 
asymmetry in the liquid and gaseous phase: fragments (liquid) appear more 
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symmetric with respect to the initial matter, while light particles (gas) are 
more neutron-rich. Hence the analyisis of the isotopic content of all reaction 
products, from pre-equilibrium emission to fragments, allows to get informa- 
tion on low-density properties of the isovector part of the nuclear interaction. 

In recent years, the properties of fragments and light clusters emitted in sys- 
tems with different initial asymmetries have been widely investigated [7,8,9,10,11,12,13], 
looking in particular at the production yields of various isotopes, as obtained 
in reactions between proton-rich and neutron-rich systems. 

More recently, the study of the isotopic content of pre-equilibrium emission has 
revealed a good sensitivity to the iso-EOS, considering the emitted neutron 
to proton ratio as a function of the kinetic energy [14]. Here we extend this 
type of investigation to fragments. Correlations between fragment charges and 
velocities have been recently observed, providing information on the interplay 
between thermal and entrance channel (collective) effects in the fragmentation 
mechanism [16,17]. Following this line, one can also investigate correlations 
between fragment isotopic content and kinematical properties, trying to get 
a deeper insight into the reaction path and to improve the understanding of 
the underlying mechanisms. In this way one can also study more in detail 
the effects of different EOS's and, in particular, of the symmetry energy on 
fragment properties [18]. 

The collision dynamics is described on the basis of approximate treatments of 
the Boltzmann-Langevin (BL) approach, i.e. of transport equations including 
fluctuations. In particular, we will adopt the Stochastic Mean Field (SMF) 
method, where fluctuations are projected onto the ordinary space [19]. 

A new method to solve the complete BL equation, by implementing the full 
fluctuations in phase space will be discussed in the last part of the manuscript 
(Section 4). 



2 The model 



Theoretically the evolution of complex systems under the influence of fluctu- 
ations can be described by a transport equation with a fluctuating term, the 
so-called Boltzmann-Langevin equation (BLE): 

f t = % + {f,H} = Icou[f]+5I[f], (1) 

where /(r, p, i) is the one-body distribution function, H(r, p, t) is the one- 
body Halmitonian and 5I[f] represents the fluctuating part of the two-body 
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collision integral [20,21]. 

Here we will follow the approximate treatment to the BLE presented in Ref. [19], 
the so-called Stochastic Mean Field (SMF) model, that consists in the imple- 
mentation of spatial density fluctuations. 

Calculations have been performed using the TWINGO code, where the test 
particle method is used to solve Eq.(l) [22]. We adopt a soft EOS, with 
compressibility modulus K = 200MeV and, for the density (p) dependence 
of the symmetry energy, we consider two representative parameterizations, 
E sym (p, I) I 'A = C sym (p)I 2 , / = (N — Z)/A : one showing a rapidly increasing 
behaviour with density, roughly proportional to p 2 (asystiff) and one with a 
kind of saturation above normal density (asysoft, SKM*) (see Ref. s [3, 23] for 
more detail). The two parameterizations obviously cross at normal density. 
The symmetry energy at densities below the normal value is larger in the 
asysoft case, while above normal density it is higher in the asystiff case. Hence 
in the low-density regime, that is the region of interest for our analysis, isospin 
effects are expected to be stronger in the asysoft case. 



3 Fragment isotopic properties and correlations 

We will focus on central collisions, b = 2 fm, considering symmetric reac- 
tions between systems having three different initial asymmetry: n2 Sn + 112 
Sn, 124 Sn + 124 Sn, 132 Sn + 132 Sn, with (N/Z) in = 1.24,1.48,1.64, respectively. 
The considered beam energy is 50 MeV/A. 1200 events have been run for each 
reaction and for each of the two symmetry energies adopted. 

The first two reactions, U2 Sn+ 112 Sn and 124 Sn+ 124 Sn have been widely inves- 
tigated both from the experimental and theoretical point of view [7,23,24,25]. 
In central collisions, after the initial collisional shock, the system expands and 
breaks up into many pieces, due to the development of volume (spinodal) and 
surface instabilities. The formation of a bubble-like configuration is observed, 
where the initial fragments are located. The average fragment multiplicity is 
approximately equal to 6 for the reactions considered here [23] . Along the re- 
action path, several nucleons are emitted at the early stage (pre-equilibrium 
emission) and/or are evaporated while fragments are formed. Primary frag- 
ments are identified by applying a coalescence procedure to the matter with 
density larger than p cut = 1/5 po (liquid phase). The remaining particles are 
considered as belonging to the gas phase. 

First, let us briefly recall some general features concerning the isotopic content 
of fragments and emitted nucleons, as obtained with the two iso-EOS's consid- 
ered. In the following we will restrict our analysis to fragments with charge in 
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Fig. 1. The N/Z of the liquid (left) and of the gas (right) phase is displayed as 
a function of the system initial N/Z. Full lines and symbols refer to the asystiff 
parameterization. Dashed lines and open symbols are for asysoft. 



the range between 3 and 10 (that we call intermediate mass fragments (IMF's) 
). The average N/Z of emitted nucleons (gas phase) and IMF's is presented 
in Fig.l as a function of the initial (N/Z) in of the three colliding systems. 
One observes that, generally, the gas phase (right panel) is more neutron- rich 
in the asysoft case, while IMF's (left panel) are more symmetric. This is due 
to the larger value of the symmetry energy at low density for the asysoft pa- 
rameterization [23]. It is interesting to note that, in the asystiff case, due to 
the rather low value of the symmetry energy, Coulomb effects dominate and 
the N/Z of the gas phase becomes lower than that for IMF's, because protons 
are preferentially emitted. Now we move to investigate in more detail correla- 
tions between fragment isotopic content and kinematical properties. The idea 
in this investigation is that fragmentation originates from the break-up of a 
composite source that expands with a given velocity field. Since neutrons and 
protons experience different forces, one may expect a different radial flow for 
the two species. In this case, the N/Z composition of the source would not be 
uniform, but would depend on the radial distance from the center or mass or, 
equivalently, on the local velocity. This trend should then be reflected in the 
fragment asymmetries. As a measure of the isotopic composition of the IMF's, 
we will consider the sum of neutrons, N = J2i^i, and protons, Z = J2i of 
all IMF's in a given kinetic energy range, in each event. Then we take the ratio 
N/Z and we consider the average over the ensemble of events. This observable 
is plotted in Fig. 2 for the three reactions. The behaviour observed is rather 
sensitive to the iso-EOS. For the proton-rich system, the N/Z decreases with 
the fragment kinetic energy, expecially in the asystiff case, where the sym- 
metry energy is relatively low at low density [23]. In this case, the Coulomb 
repulsion pushes the protons towards the surface of the system. Hence, more 
symmetric fragments acquire larger velocity. The same effects are responsible 
for the proton- rich pre-equilibrium emission observed in this case (see Fig.l). 
The decreasing trend is less pronounced in the asysoft case (right panel) be- 
cause Coulomb effects on protons are counterbalanced by the larger attraction 
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Fig. 2. The fragment N/Z (see text) as a function of the kinetic energy. Left 
panekasystiff; Rigth panel: asysoft. 

of the symmetry potential. In systems with higher initial asymmetry, the de- 
creasing trend is inversed, due to the larger neutron repulsion in neutron-rich 
systems. Larger slopes are always observed in the asysoft case. 

In conclusion, this analysis reveals the existence of significant, EOS-dependent 
correlations between the N/Z and the kinetic energy of IMF's. This correlation 
is linked to the different forces experienced by neutrons and protons along the 
fragmentation path, that in turn depend on the detail of the isovector part 
of the nuclear interaction. This study can be considered as complementary 
to pre-equilibrium emission studies [14,15]. A parallel investigation of pre- 
equilibrium and fragment emissions would be very important for a cross-check 
of model predictions against experimental observables sensitive to different 
phases of the reaction. 



4 A new way to implement the full BL equation 

The results discussed above were based on an approximate treatment of the BL 
equation, where only fluctuations in the ordinary space are considered. This is 
justified by the fact that the mult i- fragment at ion mechanism is mostly domi- 
nated by spatial density fluctuations. However, a more accurate representation 
of the full phase space dynamics would improve the description of fluctuations 
of fragment kinematical properties, allowing to investigate more sophisticated 
observables and to study the fragmentation path in deeper detail. 

Within the test partice method, that is usually adopted to solve numerically 
the BNV equation, each nucleon is represented by a collection of N test test 
particles, that are propagated according to the mean-field interaction and hard 
two-body scattering. Since collisions are treated stochastically for the single 
test particles, the dispersion around the average number of nucleon collisions 
is automatically divided by N test . Hence fluctuations are largely reduced. 
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To overcome this problem, Bauer et al. [26] proposed a method to introduce 
a correlation between close particles in phase-space. The method follows the 
idea, first applied in extended TDHF, of evolving the single-particle density 
including a perturbative mixing of Slater determinants. This is realized in the 
collision integral by means of the following procedure: The nucleon-nucleon 
(NN) cross section is reduced by a factor N test . Two test particles %\ and i 2 are 
chosen as colliding partners, and will be moved to new positions in p space 
provided that the corresponding transition probability is larger than a random 
number (extracted between and 1), as usually done in transport codes using 
the test particle method. If the two particles can collide, i.e. the final positions 
are not forbidden by the Pauli blocking, the scattering actually involves two 
"clouds" of neighbouring test particles, corresponding to two entire nucleons 
(2 x N test particles), that are chosen as the particles at closest distance in 
phase space. 

This method is conceptually simple, and moreover the implementation into 
existing transport codes is really straightforward. However, since the Pauli- 
blocking is only checked for the first two colliding test particles and not for the 
entire swarm of particles that are actually moved, at the end the fluctuation 
strength does not reproduce the expected value for fermionic systems, but 
approaches instead the classical value. 

A careful modification of the original procedure can considerably improve the 
results. The procedure that we propose can be summarized in the following 
steps: 

(i) The choice of the two colliding partners closely traces the standard recipe. 
If two test particles are allowed to collide, moving from the initial posi- 
tions pi, p 2 to new positions, p 3 , p 4 corresponding to a rotation by (9, (j>) 
with respect to the relative momentum direction, two clouds of N test par- 
ticles are moved. 

(ii) Only particles within a sphere, in coordinate space, around the center of 
mass of the two partners i\ and i 2 can belong to the clouds. The distance 
criterion is | — r c - M (i 1 ,i 2 )| < d r , where d r is a free parameter (see later 
discussion). Then a grid is introduced in momentum space, the size of 
each cell being V ce u- 

(iii) indicating by / and J the cells containing the partners i\ and i 2 , we con- 
sider the corresponding final cells I' and J' in the frame rotated by (#, 4>). 
For a given set of initial and final cells, the number of test particles that 
will be actually moved to final states is equal to the minimum between 
the occupation of the initial cells and the availability of the final ones 
n' = (1-n'): 

n t (I, I'; J, J') = min(n/, nj, n r ,nj>) 

(iv) Surrounding cells are explored with the same prescriptions, until one 
entire nucleon is reconstructed around i\ and i 2 . The search procedure is 
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symmetric with respect to the center of mass of the two partners, 
(v) A further check on the total momentum of the two clouds is performed, 
and the origin of the grid is eventually slightly displaced in order to have 
a perfect energy and momentum conservation. 

This method involves two parameters, namely the radius of the sphere in 
r space, d r , and the size of the momentum cells, V ce u. The radius d r fixes 
the spatial extension of the nucleon, that in turn influences its spreading in 
momentum space. Hence this parameter fixes somehow the extension of the 
nucleon wave packet in phase-space and generally could affect the transport 
dynamics. It can be constrained by physical arguments, and in general it 
depends on physical properties of the system, such as its dilution. The cell 
size V ce u should be small enough to allow an accurate check of the Pauli 
blocking, but large enough to contain a sufficient number of test particles to 
reduce numerical uncertainties. 

4-1 Check of fluctuations 

We discuss our procedure to build fluctuations in the context of a system of 
fermions at equilibrium, at a given density and temperature, i.e. initialized ac- 
cording to a Fermi-Dirac distribution. This is a situation that is easily reached 
in the course of a nuclear reaction, after the initial collisional shock. 

We emulate a piece of nuclear matter by taking a box with periodic boundary 
conditions. Since we will focus on fluctuations in momentum space, only one 
cell is present in coordinate space, and all particles can be chosen to collide 
(no restrictions in r space). The size of the box is I = 26 fm, and we put 2820 
nucleons, so that the density has the saturation value g = 0.16 fm~ 3 ; each 
nucleon is represented by a collection of 500 test particles. Besides, we do not 
consider any distinction between neutrons and protons, so that one nucleon 
occupies at least a phase-space volume h 3 /A. We initialize the momenta so as 
to reproduce a Fermi-Dirac profile corresponding to a temperature of 5 MeV; 
finally, we consider a constant cross section o = 160 mb. In these calculations 
the volume V ce u corresponds to a cube of size l s = 30 MeV/c (and coincides 
with V p = h 3 /(4:V r )). Calculations are stopped when the fluctuation variance 
saturates. 

First of all, we have checked that our fluctuating collision integral preserves 
the average evolution of the system, i.e does not change the one-body density 
profile. Also the average collision rate is in good agreement with analytical 
expectations [27]. 

One can directly visualize the effect of the fluctuating term by selecting a 
thin region in momentum space around the Fermi momentum. This roughly 
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corresponds to look at a bi-dimensional system. For this purpose, we adopt a 
new set of coordinates (dp, pd6, p sm(9)d(f>); for fixed p and dp, the distribution 
function will depend only on two coordinates, namely: 

f{p,8,<f>)^f{e,n) 

with di] = sin(6 l )0. We choose p = 260 MeV/c (approximately equal to the 
Fermi momentum) and Ap = 10 MeV/c. In Fig. 3 we plot the distribution 
function f(6,r]) at two different times (initial and final times). At the initial 
time, / is nearly uniform, and its fluctuations are simply due to the numerical 
noise induced by the finite number of test particles. At later times, we observe 
the growth of fluctuations, evidenced by the typical structure with "peaks and 
holes". 




Fig. 3. Distribution function f(0,rj) at time t = (left) and t = 100 fm/c (right). 



In order to check the amplitude of fluctuations built by our procedure, we 
construct the density variance in a phase- space cell of volume V, a 2 =< 
Sf(p)5f(p) >. At the Fermi surface, the expected value of the variance is 
equal to: a 2 = / eg (l — feqj/Ny = 0.25/Ny, where Ny is the number of nu- 
cleons that can be contained, at most, inside the volume V: N v = iV/h 3 . f eq 
denotes the average equilibrium value of the distribution function that, at the 
Fermi surface, is equal to 0.5. 

Since the nucleon wave packet is generally deformed, it extends over volumes 
larger than h 3 /4, hence the correct value of fluctuations is recovered in the 
"continuum limit", i.e. in large volumes V, as shown in Fig. 4. The slight 
deviation from the expected value, 0.25, is due to some difficulties, inherent 
to the adopted numerical procedure, in reconstructing entire nucleons in all 
collisional events. 
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Fig. 4. Variance of the distribution function (rescaled by Ny) as a function of the 
number Ny of nucleons contained in the cell considered to construct the fluctuation 
variance. 

5 Conclusions 

We have discussed properties of the fragmentation path in central collisions 
of charge-asymmetric systems, that can be related to the behaviour of the 
symmetry energy below normal density, thus allowing to extract information 
on fundamental quantities of the nuclear interaction. We focused on the anal- 
ysis of correlations between fragment isotopic content and kinetic energy, per- 
forming simulations based on SMF approaches. An EOS-dependent relation 
between these two observables is found. This study also allows one to get a 
deeper insight into the fragmentation mechanism. In fact, the analysis of cor- 
relations between fragment composition and velocity can be used as a clock 
of fragment formation and as an indicator of the underlying dynamics. 

A new method to improve the treatment of fluctuations in transport ap- 
proaches is also discussed. This is relevant not only for fragmentation studies, 
but in general for the dynamical description of quantum many-body systems. 
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